full_data <- readRDS('data/full_data.rds')
daily_full <- readRDS('data/daily_full.rds')

Monthly Production Data

prod_by_month <- full_data %>% 
    select(yearmonth, active_1km, monthly_oil_1km, monthly_gas_1km, 
           monthly_water_1km, monthly_boe_1km) %>%
  distinct() %>%
  filter(yearmonth <= '2023-03')

coef <- median(1032*prod_by_month$monthly_gas_1km/10^9)/median(prod_by_month$monthly_oil_1km/10^6)

full_data %>%
  ggplot(aes(x = yearmonth, group = 1)) + 
  geom_line(aes(y = monthly_oil_1km/10^6, color = 'Oil'), size=1) + 
  geom_line(aes(y = 1032*monthly_gas_1km/10^9/coef, color = 'Natural Gas'), linewidth=1) +
  scale_y_continuous(
    name = "Monthly Oil in Million Barrels",
    sec.axis = sec_axis(~.*coef, name="Monthly Gas in Billion Cubic Feet")
  ) +
  theme(
    axis.title.y = element_text(color = 'tomato'),
    axis.title.y.right = element_text(color = 'skyblue')
  ) +
  scale_color_manual(values = c("skyblue", "tomato")) +
  scale_x_discrete(breaks = prod_by_month$yearmonth[seq(1, length(prod_by_month$yearmonth), by = 7)]) +
  xlab("Date") + labs(colour="Production Type") +
    theme_minimal() +
  theme(axis.text.x = element_text(angle = 45))
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Visualizing Daily Max H2S

H2S_dm <- full_data %>%
  group_by(day, Monitor) %>%
  summarise(H2S_daily_max = max(H2S, na.rm = TRUE)) %>%
  mutate(H2S_daily_max = if_else(H2S_daily_max == -Inf, NA, H2S_daily_max))
## Warning: There were 111 warnings in `summarise()`.
## The first warning was:
## ℹ In argument: `H2S_daily_max = max(H2S, na.rm = TRUE)`.
## ℹ In group 1163: `day = 2021-01-19`, `Monitor = "G Street"`.
## Caused by warning in `max()`:
## ! no non-missing arguments to max; returning -Inf
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 110 remaining warnings.
## `summarise()` has grouped output by 'day'. You can override using the `.groups`
## argument.
H2S_ma <- H2S_dm %>%
  mutate(yearmonth = format(day, "%Y-%m")) %>%
  group_by(yearmonth, Monitor) %>%
  summarise(H2S_monthly_average = mean(H2S_daily_max, na.rm = TRUE))
## `summarise()` has grouped output by 'yearmonth'. You can override using the
## `.groups` argument.
H2S_dm_graph <- H2S_dm %>%
 ggplot(aes(x=day, y=H2S_daily_max, group=Monitor, color=Monitor)) +
   geom_line() +
   labs(title = "Daily Max H2S concentration by monitor", y = 'Daily Max H2S', x = 'time') +
   theme_minimal() +
   theme(axis.text.x = element_text(angle = 45))

ggplotly(H2S_dm_graph) %>% layout(dragmode = 'pan')
H2S_dm_graph_b <- H2S_dm %>%
 ggplot(aes(x=day, y=H2S_daily_max, group=Monitor, color=Monitor)) +
   geom_line() +
   labs(title = "Daily Max H2S concentration by monitor w/o outlier", y = 'Daily Max H2S', x = 'time') +
   scale_y_continuous(limits = c(0, 100)) +
   theme_minimal() +
   theme(axis.text.x = element_text(angle = 45))

ggplotly(H2S_dm_graph_b) %>% layout(dragmode = 'pan')
H2S_ma_graph <- H2S_ma %>%
 ggplot(aes(x=yearmonth, y=H2S_monthly_average, group=Monitor, color=Monitor)) +
   geom_line() +
   labs(title = "Monthly Average H2S concentration by monitor", y = 'Monthly Average H2S', x = 'time') +
   scale_x_discrete(breaks = unique(H2S_ma$yearmonth)[seq(1, length(unique(H2S_ma$yearmonth)), by = 6)]) +
   theme_minimal() +
   theme(axis.text.x = element_text(angle = 45))

ggplotly(H2S_ma_graph) %>% layout(dragmode = 'pan')
H2S_ma_graph_b <- H2S_ma %>%
 ggplot(aes(x=yearmonth, y=H2S_monthly_average, group=Monitor, color=Monitor)) +
   geom_line() +
   labs(title = "Monthly Average H2S concentration by monitor w/o outlier", y = 'Monthly Average H2S', x = 'time') +
   scale_x_discrete(breaks = unique(H2S_ma$yearmonth)[seq(1, length(unique(H2S_ma$yearmonth)), by = 6)]) +
  scale_y_continuous(limits=c(0, 50)) +
   theme_minimal() +
   theme(axis.text.x = element_text(angle = 45))

ggplotly(H2S_ma_graph_b) %>% layout(dragmode = 'pan')

The line with the largest deviation/peak corresponds to the 213th&Chico monitor, in 2021-10 Maybe we should start from January 2022 to start modelling to remove the extreme values, or simply removing the 213 monitor since that’s set up specifically for the event.

Visualizing Daily Average H2S

# Compute daily average
H2S_da <- full_data %>%
  group_by(day, Monitor) %>%
  summarise(H2S_daily_avg = mean(H2S, na.rm=TRUE))
## `summarise()` has grouped output by 'day'. You can override using the `.groups`
## argument.
H2S_da_graph <- H2S_da %>%
 ggplot(aes(x=day, y=H2S_daily_avg, group=Monitor, color=Monitor)) +
   geom_line() +
   labs(title = "Daily Average H2S concentration by monitor", y = 'Daily Average H2S', x = 'time') +
   theme_minimal() +
   theme(axis.text.x = element_text(angle = 45))

ggplotly(H2S_da_graph) %>% layout(dragmode = 'pan')
H2S_da_graph_b <- H2S_da %>%
 ggplot(aes(x=day, y=H2S_daily_avg, group=Monitor, color=Monitor)) +
   geom_line() +
   labs(title = "Daily Average H2S concentration by monitor w/o outlier", y = 'Daily Avearge H2S', x = 'time') +
   scale_y_continuous(limits = c(0, 100)) +
   theme_minimal() +
   theme(axis.text.x = element_text(angle = 45))

ggplotly(H2S_dm_graph_b) %>% layout(dragmode = 'pan')
full_data %>%
  polarPlot(pollutant = "H2S", col = "turbo", 
            key.position = "bottom",
            key.header = "mean H2S", 
            key.footer = NULL, title = 'hi')

Monitor/refinery Locations

#https://atlas.eia.gov/datasets/6547eda91ef84cc386e23397cf834524_22/about
st_read('shapefiles/Petroleum_Refinery.shp')
## Reading layer `Petroleum_Refinery' from data source 
##   `D:\313\Year4Uni\Reading Course\H2S-study\shapefiles\Petroleum_Refinery.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 6 features and 22 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -118.411 ymin: 33.77906 xmax: -118.2337 ymax: 33.91142
## Geodetic CRS:  WGS 84
refineries <- read_sf('shapefiles/Petroleum_Refinery.shp', layer = 'Petroleum_Refinery')
refineries <- refineries %>%
  select(Company, Site, Latitude, Longitude) %>%
  mutate(refinery = case_when(Site == 'Carson' ~ 'Marathon (Carson)',
                              Site == 'El Segundo' ~ 'Chevron El Segundo',
                              Site == 'Wilmington Asphalt Plant' ~ 'Valero Refinery',
                              Site == 'Wilmington Refinery' ~ 'Marathon (Wilmington)',
                              Site == 'Torrance' ~ 'Torrance Refinery',
                              Site == 'Wilmington' ~ 'Phillips 66 (Wilmington)'))
monitor_loc <- daily_full %>%
  select(monitor_lon, monitor_lat, Monitor) %>%
  distinct()

ref_loc <- daily_full %>%
  select(Refinery) %>% distinct()

# also load in dominguez channel
st_read('shapefiles/DominguezChannel_Carson.shp')
## Reading layer `DominguezChannel_Carson' from data source 
##   `D:\313\Year4Uni\Reading Course\H2S-study\shapefiles\DominguezChannel_Carson.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 9 features and 17 fields
## Geometry type: LINESTRING
## Dimension:     XYZ
## Bounding box:  xmin: -118.3382 ymin: 33.77701 xmax: -118.2275 ymax: 33.92881
## z_range:       zmin: 0 zmax: 0
## Geodetic CRS:  WGS 84
d_channel <- read_sf('shapefiles/DominguezChannel_Carson.shp', layer = 'DominguezChannel_Carson')

basemap <- leaflet() %>%
              addProviderTiles(providers$MtbMap) %>%
              addProviderTiles(providers$Stamen.TerrainBackground,
                               options = providerTileOptions(opacity = 0.9)) %>%
              addProviderTiles(providers$Stamen.TonerLines,
                options = providerTileOptions(opacity = 0.25)) %>%
              addProviderTiles(providers$Stamen.TonerLabels) %>%
              addCircleMarkers(data = monitor_loc, 
                         lng=~monitor_lon, lat=~monitor_lat, 
                         fillOpacity = 0.85, radius = 4,
                         fillColor = , weight = 1, stroke = T, col = '#BA55D3',
                         label = ~Monitor) %>%
              addCircleMarkers(data = refineries, 
                         lng=~Longitude, lat=~Latitude, 
                         fillOpacity = 0.95, radius = 6,
                         fillColor = , weight = 1, stroke = T, col = '#8B0000',
                         label = ~refinery)


basemap
saveWidget(basemap, file="basemap.html")

Since Feb 2022

# Create a polar map
# TBD: include markers for refineries
polarmap1 <- polarMap(
  daily_full[complete.cases(daily_full),] %>%
    filter(day >= '2022-02-01',
           !(Monitor %in% c('West HS', 'Elm Ave', 'North HS'))) %>%
    rename(wd = wd_avg,
           ws = ws_avg),
  pollutant = "H2S_daily_avg",
  latitude = "monitor_lat",
  longitude = "monitor_lon",
  popup = "Monitor",
  provider = "Stamen.Toner",
  d.icon = 150,
  d.fig = 2.5,
  alpha = 0.95,
  key = TRUE
)

polarmap2 <- polarMap(
  daily_full[complete.cases(daily_full),] %>%
    filter(day >= '2022-02-01',
           Monitor %in% c('West HS', 'Elm Ave', 'North HS')) %>%
    rename(wd = wd_avg,
           ws = ws_avg),
  pollutant = "H2S_daily_avg",
  latitude = "monitor_lat",
  longitude = "monitor_lon",
  popup = "Monitor",
  provider = "Stamen.Toner",
  d.icon = 150,
  d.fig = 2.5,
  alpha = 0.95,
  key = TRUE,
  statistic = 'nwr'
)
## 
Creating Polar Markers â– â– â– â– â– â– â– â– â– â– â–  33% | ETA: 40s

Creating Polar Markers
## â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â–  67% | ETA: 21s


polarmap1
polarmap2
sincefeb2022_mon <- daily_full[complete.cases(daily_full),] %>%
    filter(day >= '2022-02-01') %>% pull(Monitor) %>% unique()

sincefeb2022_data <- daily_full[complete.cases(daily_full),] %>%
    filter(day >= '2022-02-01') %>%
    rename(wd = wd_avg,
           ws = ws_avg)

for (mon in sincefeb2022_mon) {
  polarPlot(sincefeb2022_data %>% filter(Monitor == mon),
    pollutant = "H2S_daily_avg",
    type = 'Monitor',
    latitude = "monitor_lat",
    longitude = "monitor_lon",
    k = 30)
}

polarPlot(daily_full[complete.cases(daily_full),] %>%
    filter(day >= '2022-02-01',
           Monitor == 'West HS') %>%
    rename(wd = wd_avg,
           ws = ws_avg),
    pollutant = "H2S_daily_avg",
    latitude = "monitor_lat",
    longitude = "monitor_lon",
    k = 10)

polarPlot(daily_full[complete.cases(daily_full),] %>%
    filter(day >= '2022-02-01',
           Monitor == 'North HS') %>%
    rename(wd = wd_avg,
           ws = ws_avg),
    pollutant = "H2S_daily_avg",
    latitude = "monitor_lat",
    longitude = "monitor_lon",
    k = 10)

polarPlot(daily_full[complete.cases(daily_full),] %>%
    filter(day >= '2022-02-01',
           Monitor == 'Elm Ave') %>%
    rename(wd = wd_avg,
           ws = ws_avg),
    pollutant = "H2S_daily_avg",
    latitude = "monitor_lat",
    longitude = "monitor_lon",
    k = 10)

Disaster Period

dis_mon <- daily_full[complete.cases(daily_full),] %>%
    filter(year == '2021', month %in% c('10', '11', '12')) %>% pull(Monitor) %>% unique()

dis_data <- daily_full[complete.cases(daily_full),] %>%
    filter(year == '2021', month %in% c('10', '11', '12')) %>%
    rename(wd = wd_avg,
           ws = ws_avg)

for (mon in dis_mon) {
  polarPlot(dis_data %>% filter(Monitor == mon),
    pollutant = "H2S_daily_avg",
    type = 'Monitor',
    latitude = "monitor_lat",
    longitude = "monitor_lon",
    k = 30)
}

# save these plots
getplot <- function(mon) {
  mon_plot <- polarPlot(dis_data %>% filter(Monitor == mon),
    pollutant = "H2S_daily_avg",
    latitude = "monitor_lat",
    longitude = "monitor_lon",
    k = 30,
    key = FALSE)
  
  return(mon_plot)
}

# for loop doesn't work with png() and bg='transparent' for some reason...
mon <- dis_mon[1]
mon_plot <- getplot(mon)

mon <- str_replace_all(str_replace_all(mon, '[()]', ''), ' ', '_')
png(paste0("polarplots/disaster/", mon, ".png"), width=20, height=15, units="cm", res=800, bg="transparent")
mon_plot
## 
## ── openair object ──────────────────────────────────────────────────────────────
## Created by:
## â–¶ `polarPlot(mydata = dis_data %>% filter(Monitor == mon), pollutant =
## "H2S_daily_avg", ` and ` k = 30, key = FALSE, latitude = "monitor_lat",
## longitude = "monitor_lon")`
## 
## 
## ── This contains: ──
## 
## 
## 
## A single data frame:
## • $data [with no subset structure]
## 
## 
## 
## A single plot:
## • $plot [with no subset structure]
dev.off()
## png 
##   2
mon <- dis_mon[2]
mon_plot <- getplot(mon)

mon <- str_replace_all(str_replace_all(mon, '[()]', ''), ' ', '_')
png(paste0("polarplots/disaster/", mon, ".png"), width=20, height=15, units="cm", res=800, bg="transparent")
mon_plot
## ── openair object ──────────────────────────────────────────────────────────────
## Created by:â–¶ `polarPlot(mydata = dis_data %>% filter(Monitor == mon), pollutant =
## "H2S_daily_avg", ` and ` k = 30, key = FALSE, latitude = "monitor_lat",
## longitude = "monitor_lon")`
## ── This contains: ──
## 
## A single data frame:• $data [with no subset structure]
## 
## A single plot:• $plot [with no subset structure]
dev.off()
## png 
##   2
mon <- dis_mon[3]
mon_plot <- getplot(mon)

mon <- str_replace_all(str_replace_all(mon, '[()]', ''), ' ', '_')
png(paste0("polarplots/disaster/", mon, ".png"), width=20, height=15, units="cm", res=800, bg="transparent")
mon_plot
## ── openair object ──────────────────────────────────────────────────────────────
## Created by:â–¶ `polarPlot(mydata = dis_data %>% filter(Monitor == mon), pollutant =
## "H2S_daily_avg", ` and ` k = 30, key = FALSE, latitude = "monitor_lat",
## longitude = "monitor_lon")`
## ── This contains: ──
## 
## A single data frame:• $data [with no subset structure]
## 
## A single plot:• $plot [with no subset structure]
dev.off()
## png 
##   2
mon <- dis_mon[4]
mon_plot <- getplot(mon)

mon <- str_replace_all(str_replace_all(mon, '[()]', ''), ' ', '_')
png(paste0("polarplots/disaster/", mon, ".png"), width=20, height=15, units="cm", res=800, bg="transparent")
mon_plot
## ── openair object ──────────────────────────────────────────────────────────────
## Created by:â–¶ `polarPlot(mydata = dis_data %>% filter(Monitor == mon), pollutant =
## "H2S_daily_avg", ` and ` k = 30, key = FALSE, latitude = "monitor_lat",
## longitude = "monitor_lon")`
## ── This contains: ──
## 
## A single data frame:• $data [with no subset structure]
## 
## A single plot:• $plot [with no subset structure]
dev.off()
## png 
##   2
mon <- dis_mon[5]
mon_plot <- getplot(mon)

mon <- str_replace_all(str_replace_all(mon, '[()]', ''), ' ', '_')
png(paste0("polarplots/disaster/", mon, ".png"), width=20, height=15, units="cm", res=800, bg="transparent")
mon_plot
## ── openair object ──────────────────────────────────────────────────────────────
## Created by:â–¶ `polarPlot(mydata = dis_data %>% filter(Monitor == mon), pollutant =
## "H2S_daily_avg", ` and ` k = 30, key = FALSE, latitude = "monitor_lat",
## longitude = "monitor_lon")`
## ── This contains: ──
## 
## A single data frame:• $data [with no subset structure]
## 
## A single plot:• $plot [with no subset structure]
dev.off()
## png 
##   2
mon <- dis_mon[6]
mon_plot <- getplot(mon)

mon <- str_replace_all(str_replace_all(mon, '[()]', ''), ' ', '_')
png(paste0("polarplots/disaster/", mon, ".png"), width=20, height=15, units="cm", res=800, bg="transparent")
mon_plot
## ── openair object ──────────────────────────────────────────────────────────────
## Created by:â–¶ `polarPlot(mydata = dis_data %>% filter(Monitor == mon), pollutant =
## "H2S_daily_avg", ` and ` k = 30, key = FALSE, latitude = "monitor_lat",
## longitude = "monitor_lon")`
## ── This contains: ──
## 
## A single data frame:• $data [with no subset structure]
## 
## A single plot:• $plot [with no subset structure]
dev.off()
## png 
##   2
mon <- dis_mon[7]
mon_plot <- getplot(mon)

mon <- str_replace_all(str_replace_all(mon, '[()]', ''), ' ', '_')
png(paste0("polarplots/disaster/", mon, ".png"), width=20, height=15, units="cm", res=800, bg="transparent")
mon_plot
## ── openair object ──────────────────────────────────────────────────────────────
## Created by:â–¶ `polarPlot(mydata = dis_data %>% filter(Monitor == mon), pollutant =
## "H2S_daily_avg", ` and ` k = 30, key = FALSE, latitude = "monitor_lat",
## longitude = "monitor_lon")`
## ── This contains: ──
## 
## A single data frame:• $data [with no subset structure]
## 
## A single plot:• $plot [with no subset structure]
dev.off()
## png 
##   2
mon <- dis_mon[8]
mon_plot <- getplot(mon)

mon <- str_replace_all(str_replace_all(mon, '[()]', ''), ' ', '_')
png(paste0("polarplots/disaster/", mon, ".png"), width=20, height=15, units="cm", res=800, bg="transparent")
mon_plot
## ── openair object ──────────────────────────────────────────────────────────────
## Created by:â–¶ `polarPlot(mydata = dis_data %>% filter(Monitor == mon), pollutant =
## "H2S_daily_avg", ` and ` k = 30, key = FALSE, latitude = "monitor_lat",
## longitude = "monitor_lon")`
## ── This contains: ──
## 
## A single data frame:• $data [with no subset structure]
## 
## A single plot:• $plot [with no subset structure]
dev.off()
## png 
##   2
mon <- dis_mon[9]
mon_plot <- getplot(mon)

mon <- str_replace_all(str_replace_all(mon, '[()]', ''), ' ', '_')
png(paste0("polarplots/disaster/", mon, ".png"), width=20, height=15, units="cm", res=800, bg="transparent")
mon_plot
## ── openair object ──────────────────────────────────────────────────────────────
## Created by:â–¶ `polarPlot(mydata = dis_data %>% filter(Monitor == mon), pollutant =
## "H2S_daily_avg", ` and ` k = 30, key = FALSE, latitude = "monitor_lat",
## longitude = "monitor_lon")`
## ── This contains: ──
## 
## A single data frame:• $data [with no subset structure]
## 
## A single plot:• $plot [with no subset structure]
dev.off()
## png 
##   2
mon <- dis_mon[10]
mon_plot <- getplot(mon)

mon <- str_replace_all(str_replace_all(mon, '[()]', ''), ' ', '_')
png(paste0("polarplots/disaster/", mon, ".png"), width=20, height=15, units="cm", res=800, bg="transparent")
mon_plot
## ── openair object ──────────────────────────────────────────────────────────────
## Created by:â–¶ `polarPlot(mydata = dis_data %>% filter(Monitor == mon), pollutant =
## "H2S_daily_avg", ` and ` k = 30, key = FALSE, latitude = "monitor_lat",
## longitude = "monitor_lon")`
## ── This contains: ──
## 
## A single data frame:• $data [with no subset structure]
## 
## A single plot:• $plot [with no subset structure]
dev.off()
## png 
##   2
mon <- dis_mon[11]
mon_plot <- getplot(mon)

mon <- str_replace_all(str_replace_all(mon, '[()]', ''), ' ', '_')
png(paste0("polarplots/disaster/", mon, ".png"), width=20, height=15, units="cm", res=800, bg="transparent")
mon_plot
## ── openair object ──────────────────────────────────────────────────────────────
## Created by:â–¶ `polarPlot(mydata = dis_data %>% filter(Monitor == mon), pollutant =
## "H2S_daily_avg", ` and ` k = 30, key = FALSE, latitude = "monitor_lat",
## longitude = "monitor_lon")`
## ── This contains: ──
## 
## A single data frame:• $data [with no subset structure]
## 
## A single plot:• $plot [with no subset structure]
dev.off()
## png 
##   2
mon <- dis_mon[12]
mon_plot <- getplot(mon)

mon <- str_replace_all(str_replace_all(mon, '[()]', ''), ' ', '_')
png(paste0("polarplots/disaster/", mon, ".png"), width=20, height=15, units="cm", res=800, bg="transparent")
mon_plot
## ── openair object ──────────────────────────────────────────────────────────────
## Created by:â–¶ `polarPlot(mydata = dis_data %>% filter(Monitor == mon), pollutant =
## "H2S_daily_avg", ` and ` k = 30, key = FALSE, latitude = "monitor_lat",
## longitude = "monitor_lon")`
## ── This contains: ──
## 
## A single data frame:• $data [with no subset structure]
## 
## A single plot:• $plot [with no subset structure]
dev.off()
## png 
##   2
mon <- dis_mon[13]
mon_plot <- getplot(mon)

mon <- str_replace_all(str_replace_all(mon, '[()]', ''), ' ', '_')
png(paste0("polarplots/disaster/", mon, ".png"), width=20, height=15, units="cm", res=800, bg="transparent")
mon_plot
## ── openair object ──────────────────────────────────────────────────────────────
## Created by:â–¶ `polarPlot(mydata = dis_data %>% filter(Monitor == mon), pollutant =
## "H2S_daily_avg", ` and ` k = 30, key = FALSE, latitude = "monitor_lat",
## longitude = "monitor_lon")`
## ── This contains: ──
## 
## A single data frame:• $data [with no subset structure]
## 
## A single plot:• $plot [with no subset structure]
dev.off()
## png 
##   2
parOrig <- par()
par(parOrig)
## Warning in par(parOrig): graphical parameter "cin" cannot be set
## Warning in par(parOrig): graphical parameter "cra" cannot be set
## Warning in par(parOrig): graphical parameter "csi" cannot be set
## Warning in par(parOrig): graphical parameter "cxy" cannot be set
## Warning in par(parOrig): graphical parameter "din" cannot be set
## Warning in par(parOrig): graphical parameter "page" cannot be set
plot (1:3)

Since 2020

Full

polarmap3 <- polarMap(
  daily_full[complete.cases(daily_full),] %>%
    rename(wd = wd_avg,
           ws = ws_avg),
  pollutant = "H2S_daily_avg",
  latitude = "monitor_lat",
  longitude = "monitor_lon",
  popup = "Monitor",
  provider = "Stamen.Toner",
  d.icon = 150,
  d.fig = 2.5,
  alpha = 0.95,
  key = TRUE
)
## 
Creating Polar Markers â– â– â– â– â– â– â– â– â– â– â– â– â– â– â–  46% | ETA: 1s

Creating Polar Markers
## â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â–  69% | ETA: 1s

Creating Polar Markers
## â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â–  77% | ETA: 1s

Creating Polar Markers
## â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â–  85% | ETA: 1s

Creating Polar Markers
## â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â–  92% | ETA: 0s


## Warning: There was 1 warning in `dplyr::mutate()`.
## ℹ In argument: `plot = purrr::map(data, fun, .progress = "Creating Polar
##   Markers")`.
## Caused by warning:
## ! Not enough data to fit surface.
## Try reducing the value of the smoothing parameter, k to less than 100. 
## Or use statistic = 'nwr'.
polarmap3
polarPlot(daily_full[complete.cases(daily_full),] %>%
    rename(wd = wd_avg,
           ws = ws_avg) %>% 
      filter(Monitor == '213th & Chico'),
    pollutant = "H2S_daily_avg",
    type = 'Monitor',
    latitude = "monitor_lat",
    longitude = "monitor_lon",
    k = 30)

full_mon <- daily_full[complete.cases(daily_full),] %>% 
  pull(Monitor) %>% unique()

full_data <- daily_full[complete.cases(daily_full),] %>%
    rename(wd = wd_avg,
           ws = ws_avg)

for (mon in full_mon) {
  polarPlot(full_data %>% filter(Monitor == mon),
    pollutant = "H2S_daily_avg",
    type = 'Monitor',
    latitude = "monitor_lat",
    longitude = "monitor_lon",
    k = 30)
}

Normal period

polarmap4 <- polarMap(
  daily_full[complete.cases(daily_full),] %>%
    filter(!(year == '2021' & month %in% c('10', '11', '12'))) %>%
    rename(wd = wd_avg,
           ws = ws_avg),
  pollutant = "H2S_daily_avg",
  latitude = "monitor_lat",
  longitude = "monitor_lon",
  popup = "Monitor",
  provider = "Stamen.Toner",
  d.icon = 150,
  d.fig = 2.5,
  alpha = 0.95,
  key = TRUE
)
## 
Creating Polar Markers â– â– â– â– â– â– â– â– â– â– â– â– â– â– â–  46% | ETA: 1s

Creating Polar Markers
## â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â–  69% | ETA: 1s

Creating Polar Markers
## â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â–  77% | ETA: 1s

Creating Polar Markers
## â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â–  85% | ETA: 0s

Creating Polar Markers
## â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â– â–  92% | ETA: 0s


## Warning: There was 1 warning in `dplyr::mutate()`.
## ℹ In argument: `plot = purrr::map(data, fun, .progress = "Creating Polar
##   Markers")`.
## Caused by warning:
## ! Not enough data to fit surface.
## Try reducing the value of the smoothing parameter, k to less than 100. 
## Or use statistic = 'nwr'.
polarmap4
norm_mon <- daily_full[complete.cases(daily_full),] %>% 
  filter(!(year == '2021' & month %in% c('10', '11', '12'))) %>%
  pull(Monitor) %>% unique()

norm_data <- daily_full[complete.cases(daily_full),] %>%
  filter(!(year == '2021' & month %in% c('10', '11', '12'))) %>%
  rename(wd = wd_avg,
           ws = ws_avg)

# 213 needs lower k due to lack of data
polar_norm_213 <- polarPlot(norm_data %>% filter(Monitor == '213th & Chico'),
  pollutant = "H2S_daily_avg",
  latitude = "monitor_lat",
  longitude = "monitor_lon",
  k = 10,
  key = FALSE)

png("polarplots/normal/213th_&_Chico.png", width=20, height=15, units="cm", res=800, bg="transparent")
polar_norm_213
## 
## ── openair object ──────────────────────────────────────────────────────────────
## Created by:
## â–¶ `polarPlot(mydata = norm_data %>% filter(Monitor == "213th & Chico"), `, `
## pollutant = "H2S_daily_avg", k = 10, key = FALSE, latitude = "monitor_lat", `,
## and ` longitude = "monitor_lon")`
## 
## 
## ── This contains: ──
## 
## 
## 
## A single data frame:
## • $data [with no subset structure]
## 
## 
## 
## A single plot:
## • $plot [with no subset structure]
dev.off()
## png 
##   2
plot(polar_norm_213, bg = 'transparent')

dev.copy(png,'myplot.png',  bg = 'transparent')
## png 
##   3
# for the rest, use k=30
getplot <- function(mon) {
  mon_plot <- polarPlot(norm_data %>% filter(Monitor == mon),
    pollutant = "H2S_daily_avg",
    latitude = "monitor_lat",
    longitude = "monitor_lon",
    #limits = c(0, 1.5),
    k = 30,
    key = FALSE)
  
  return(mon_plot)
}

# for loop doesn't work with png() and bg='transparent' for some reason...
mon <- norm_mon[2]
mon_plot <- getplot(mon)
mon <- str_replace_all(str_replace_all(mon, '[()]', ''), ' ', '_')
png(paste0("polarplots/normal/", mon, ".png"), width=20, height=15, units="cm", res=800, bg="transparent")
mon_plot
## ── openair object ──────────────────────────────────────────────────────────────
## Created by:â–¶ `polarPlot(mydata = norm_data %>% filter(Monitor == mon), pollutant =
## "H2S_daily_avg", ` and ` k = 30, key = FALSE, latitude = "monitor_lat",
## longitude = "monitor_lon")`
## ── This contains: ──
## 
## A single data frame:• $data [with no subset structure]
## 
## A single plot:• $plot [with no subset structure]
dev.off()
## png 
##   2
mon <- norm_mon[3]
mon_plot <- getplot(mon)

mon <- str_replace_all(str_replace_all(mon, '[()]', ''), ' ', '_')
png(paste0("polarplots/normal/", mon, ".png"), width=20, height=15, units="cm", res=800, bg="transparent")
mon_plot
## ── openair object ──────────────────────────────────────────────────────────────
## Created by:â–¶ `polarPlot(mydata = norm_data %>% filter(Monitor == mon), pollutant =
## "H2S_daily_avg", ` and ` k = 30, key = FALSE, latitude = "monitor_lat",
## longitude = "monitor_lon")`
## ── This contains: ──
## 
## A single data frame:• $data [with no subset structure]
## 
## A single plot:• $plot [with no subset structure]
dev.off()
## png 
##   2
mon <- norm_mon[4]
mon_plot <- getplot(mon)

mon <- str_replace_all(str_replace_all(mon, '[()]', ''), ' ', '_')
png(paste0("polarplots/normal/", mon, ".png"), width=20, height=15, units="cm", res=800, bg="transparent")
mon_plot
## ── openair object ──────────────────────────────────────────────────────────────
## Created by:â–¶ `polarPlot(mydata = norm_data %>% filter(Monitor == mon), pollutant =
## "H2S_daily_avg", ` and ` k = 30, key = FALSE, latitude = "monitor_lat",
## longitude = "monitor_lon")`
## ── This contains: ──
## 
## A single data frame:• $data [with no subset structure]
## 
## A single plot:• $plot [with no subset structure]
dev.off()
## png 
##   2
mon <- norm_mon[5]
mon_plot <- getplot(mon)

mon <- str_replace_all(str_replace_all(mon, '[()]', ''), ' ', '_')
png(paste0("polarplots/normal/", mon, ".png"), width=20, height=15, units="cm", res=800, bg="transparent")
mon_plot
## ── openair object ──────────────────────────────────────────────────────────────
## Created by:â–¶ `polarPlot(mydata = norm_data %>% filter(Monitor == mon), pollutant =
## "H2S_daily_avg", ` and ` k = 30, key = FALSE, latitude = "monitor_lat",
## longitude = "monitor_lon")`
## ── This contains: ──
## 
## A single data frame:• $data [with no subset structure]
## 
## A single plot:• $plot [with no subset structure]
dev.off()
## png 
##   2
mon <- norm_mon[6]
mon_plot <- getplot(mon)

mon <- str_replace_all(str_replace_all(mon, '[()]', ''), ' ', '_')
png(paste0("polarplots/normal/", mon, ".png"), width=20, height=15, units="cm", res=800, bg="transparent")
mon_plot
## ── openair object ──────────────────────────────────────────────────────────────
## Created by:â–¶ `polarPlot(mydata = norm_data %>% filter(Monitor == mon), pollutant =
## "H2S_daily_avg", ` and ` k = 30, key = FALSE, latitude = "monitor_lat",
## longitude = "monitor_lon")`
## ── This contains: ──
## 
## A single data frame:• $data [with no subset structure]
## 
## A single plot:• $plot [with no subset structure]
dev.off()
## png 
##   2
mon <- norm_mon[7]
mon_plot <- getplot(mon)

mon <- str_replace_all(str_replace_all(mon, '[()]', ''), ' ', '_')
png(paste0("polarplots/normal/", mon, ".png"), width=20, height=15, units="cm", res=800, bg="transparent")
mon_plot
## ── openair object ──────────────────────────────────────────────────────────────
## Created by:â–¶ `polarPlot(mydata = norm_data %>% filter(Monitor == mon), pollutant =
## "H2S_daily_avg", ` and ` k = 30, key = FALSE, latitude = "monitor_lat",
## longitude = "monitor_lon")`
## ── This contains: ──
## 
## A single data frame:• $data [with no subset structure]
## 
## A single plot:• $plot [with no subset structure]
dev.off()
## png 
##   2
mon <- norm_mon[8]
mon_plot <- getplot(mon)

mon <- str_replace_all(str_replace_all(mon, '[()]', ''), ' ', '_')
png(paste0("polarplots/normal/", mon, ".png"), width=20, height=15, units="cm", res=800, bg="transparent")
mon_plot
## ── openair object ──────────────────────────────────────────────────────────────
## Created by:â–¶ `polarPlot(mydata = norm_data %>% filter(Monitor == mon), pollutant =
## "H2S_daily_avg", ` and ` k = 30, key = FALSE, latitude = "monitor_lat",
## longitude = "monitor_lon")`
## ── This contains: ──
## 
## A single data frame:• $data [with no subset structure]
## 
## A single plot:• $plot [with no subset structure]
dev.off()
## png 
##   2
mon <- norm_mon[9]
mon_plot <- getplot(mon)

mon <- str_replace_all(str_replace_all(mon, '[()]', ''), ' ', '_')
png(paste0("polarplots/normal/", mon, ".png"), width=20, height=15, units="cm", res=800, bg="transparent")
mon_plot
## ── openair object ──────────────────────────────────────────────────────────────
## Created by:â–¶ `polarPlot(mydata = norm_data %>% filter(Monitor == mon), pollutant =
## "H2S_daily_avg", ` and ` k = 30, key = FALSE, latitude = "monitor_lat",
## longitude = "monitor_lon")`
## ── This contains: ──
## 
## A single data frame:• $data [with no subset structure]
## 
## A single plot:• $plot [with no subset structure]
dev.off()
## png 
##   2
mon <- norm_mon[10]
mon_plot <- getplot(mon)

mon <- str_replace_all(str_replace_all(mon, '[()]', ''), ' ', '_')
png(paste0("polarplots/normal/", mon, ".png"), width=20, height=15, units="cm", res=800, bg="transparent")
mon_plot
## ── openair object ──────────────────────────────────────────────────────────────
## Created by:â–¶ `polarPlot(mydata = norm_data %>% filter(Monitor == mon), pollutant =
## "H2S_daily_avg", ` and ` k = 30, key = FALSE, latitude = "monitor_lat",
## longitude = "monitor_lon")`
## ── This contains: ──
## 
## A single data frame:• $data [with no subset structure]
## 
## A single plot:• $plot [with no subset structure]
dev.off()
## png 
##   2
mon <- norm_mon[11]
mon_plot <- getplot(mon)

mon <- str_replace_all(str_replace_all(mon, '[()]', ''), ' ', '_')
png(paste0("polarplots/normal/", mon, ".png"), width=20, height=15, units="cm", res=800, bg="transparent")
mon_plot
## ── openair object ──────────────────────────────────────────────────────────────
## Created by:â–¶ `polarPlot(mydata = norm_data %>% filter(Monitor == mon), pollutant =
## "H2S_daily_avg", ` and ` k = 30, key = FALSE, latitude = "monitor_lat",
## longitude = "monitor_lon")`
## ── This contains: ──
## 
## A single data frame:• $data [with no subset structure]
## 
## A single plot:• $plot [with no subset structure]
dev.off()
## png 
##   2
mon <- norm_mon[12]
mon_plot <- getplot(mon)

mon <- str_replace_all(str_replace_all(mon, '[()]', ''), ' ', '_')
png(paste0("polarplots/normal/", mon, ".png"), width=20, height=15, units="cm", res=800, bg="transparent")
mon_plot
## ── openair object ──────────────────────────────────────────────────────────────
## Created by:â–¶ `polarPlot(mydata = norm_data %>% filter(Monitor == mon), pollutant =
## "H2S_daily_avg", ` and ` k = 30, key = FALSE, latitude = "monitor_lat",
## longitude = "monitor_lon")`
## ── This contains: ──
## 
## A single data frame:• $data [with no subset structure]
## 
## A single plot:• $plot [with no subset structure]
dev.off()
## png 
##   2
mon <- norm_mon[13]
mon_plot <- getplot(mon)

mon <- str_replace_all(str_replace_all(mon, '[()]', ''), ' ', '_')
png(paste0("polarplots/normal/", mon, ".png"), width=20, height=15, units="cm", res=800, bg="transparent")
mon_plot
## ── openair object ──────────────────────────────────────────────────────────────
## Created by:â–¶ `polarPlot(mydata = norm_data %>% filter(Monitor == mon), pollutant =
## "H2S_daily_avg", ` and ` k = 30, key = FALSE, latitude = "monitor_lat",
## longitude = "monitor_lon")`
## ── This contains: ──
## 
## A single data frame:• $data [with no subset structure]
## 
## A single plot:• $plot [with no subset structure]
dev.off()
## png 
##   2
parOrig <- par()
par(parOrig)
## Warning in par(parOrig): graphical parameter "cin" cannot be set
## Warning in par(parOrig): graphical parameter "cra" cannot be set
## Warning in par(parOrig): graphical parameter "csi" cannot be set
## Warning in par(parOrig): graphical parameter "cxy" cannot be set
## Warning in par(parOrig): graphical parameter "din" cannot be set
## Warning in par(parOrig): graphical parameter "page" cannot be set
plot (1:3)

Visualize Odor Complaint data

odor <- daily_full %>%
  select(day, county, num_odor_complaints) %>%
  distinct()

odor_graph <- odor %>%
   ggplot(aes(x=day, y=num_odor_complaints, group=county, color=county)) +
   geom_line() +
   labs(title = "Number of odor complaints over time", y = 'Number of odor complaints', x = 'time') +
   theme_minimal() +
   theme(axis.text.x = element_text(angle = 45))

ggplotly(odor_graph) %>% layout(dragmode = 'pan')
odor_graph_b <- odor %>%
   ggplot(aes(x=day, y=num_odor_complaints, group=county, color=county)) +
   geom_line() +
   labs(title = "Number of odor complaints over time", y = 'Number of odor complaints', x = 'time') +
  scale_y_continuous(limits = c(0, 30)) +
   theme_minimal() +
   theme(axis.text.x = element_text(angle = 45))

ggplotly(odor_graph_b) %>% layout(dragmode = 'pan')